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Abstract 

We consider a simple model consisting of particles with four bonding sites ("patches"), two of 
type A and two of type B, on the square lattice, and investigate its global phase behavior by 
simulations and theory. We set the interaction between B patches to zero and calculate the phase 
diagram as the ratio between the AB and the AA interactions, e^^, varies. In line with previous 
work, on three-dimensional off-lattice models, we show that the liquid-vapor phase diagram exhibits 
a re-entrant or "pinched" shape for the same range of e^^, suggesting that the ratio of the energy 
scales - and the corresponding empty fluid regime - is independent of the dimensionality of the 
system and of the lattice structure. In addition, the model exhibits an order-disorder transition 
that is ferromagnetic in the re-entrant regime. The use of low-dimensional lattice models allows 
the simulation of sufficiently large systems to establish the nature of the liquid- vapor critical points 
and to describe the structure of the liquid phase in the empty fluid regime, where the size of the 
"voids" increases as the temperature decreases. We have found that the liquid-vapor critical point 
is in the 2D Ising universality class, with a scaling region that decreases rapidly as the temperature 
decreases. 

The results of simulations and theoretical analysis suggest that the line of order-disorder transi- 
tions intersects the condensation line at a multicritical point at zero temperature and density, for 

patchy particle models with a re-entrant, empty fluid, regime. 

PACS numbers: 64.60Cn, 61.20.Gy 



2 



I. INTRODUCTION 



One of the scientific and technological revolutions currently in progress is the increas- 
ing ability to miniaturize material design and manufacturing components. Advances in 
the chemical synthesis and fabrication of nanometer-to-micrometer sized particles have pro- 
duced a variety of new particles. Their organization into more complex structures remains, 
however, a great challenge. A promising approach inspired by Nature is nanoparticle self- 
assembly. The structure of the self-assembled clusters, which range from chains to rings and 
complex branched structures, depends crucially on the anisotropy of the particle shapes and 
interactions and may compete with the clustering that drives condensation, giving rise to 



novel macroscopic behavior 



Indeed, patchy particle models with dissimilar patches {A and B) were recently intro- 
duced in this context and revealed that the criticality depends on the type of clusters that 



are formed 



Of particular interest are systems where the self-assembled clusters are 



long linear chains connected by junctions, as the liquid- vapor transition of these network 
(percolated) fluids may be viewed as the condensation of these junctions In addition 

to ferrofluids or electro-rheological fluids, colloids with distinct patchy interactions may be 



synthesized by the selective functionalisation of specific areas of the particles 

Primitive models of patchy particles with identical js, 9| and distinct patches share the 
physics of limited valence materials, namely the existence of stable liquid states of van- 
ishingly small density (empty liquids), and provide a route to equilibrium gels lOj. In 
addition, distinct-patch models allow a unique control of the effective valence through the 
temperature T. In three-dimensional (3D) off-lattice models consisting of particles with two 
types of patches, A and B, where the interaction between B patches is set to zero, the 
topology of the liquid-vapor diagram is determined by the ratio between the AB and the 
A A interactions, e\^. As e^^ decreases in the range | < < |, the low-temperature 
liquid- vapor coexistence region also decreases The binodal exhibits a characteristic 



re-entrant or "pinched" s iape with the coexisting liquid density vanishing as the tempera- 



ture approaches zero [11 



^AB ~ I there is no condensation, and above e^^ 



there is no re-entrant behavior [13[. Both the scaling of the vanishing critical parameters 
and the re-entrant phase behavior are predicted correctly by Wertheim's thermodynamic 
first-order perturbation theory [ill, Q, [l4, 15|. The theory also reveals that the re-entrant 
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phase behavior is driven by the balance of two entropic contributions: the higher entropy of 
the junctions and the lower entropy of the chains in the (network) liquid phase, as suggested 
a decade ago on the basis of a hierarchical theory of network fluids [16| . 

The feature that makes patchy particles ideally suited to the investigation of the interplay 
between self-assembly and condensation is the fact that both the thermodynamic and struc- 
tural properties of patchy particle systems can be predicted with a high degree of accuracy 
by the thermodynamic perturbation theory of Wertheim and the Flory-Stockmayer theory 



of polymerization 17H19|. It is then possible to study the phase behavior of patchy parti- 



cles using reliable liquid-state theories and to use this knowledge to design the models and 
guide the simulations, the results of which validated the theoretical predictions js, 9, 11, 12 1. 



There remain, however, two open questions: 1. What is the nature of the liquid- vapor 
critical point, in models with an empty fluid regime ? and 2. Are there ordered phases 
that pre-empt the empty fluid regime or. What is the topology of the global phase diagram 
? These are difficult questions that will be addressed here by considering simple patchy 
particle models on the square lattice. 

In systems with two bonding sites per particle, only (polydisperse) linear chains form 
and there is no liquid- vapor phase transition 2oi. If the chains are sufficiently stiff they 
undergo an ordering transition at fixed concentration, as the temperature decreases below 
the bonding temperature. The interplay between the self-assembly process, driven by the 
bonding interactions, and the ordering transition, driven by the anisotropic shape of the 
bonded clusters has been investigated for a two-dimensional (2D) model consisting of par- 
ticles with two bonding sites, on the square lattice (self-assembling rigid rods or SARR 
model). It was shown that bonding drives ordering and that the ordering enhances bonding 



2l| . Subsequently, extensive Monte Carlo simulations were carried out to investigate the 



nature of the ordering transition that was shown to be in the Ising 2D universality class, as 
in models where the rods are mono disperse 22|. The scaling region, however, was found to 
depend strongly on the temperature |22l. |23|. 

In this paper, we consider the 2A2B model consisting of particles with four patches, two of 
type A and two of type B, on the square lattice and investigate its global phase behavior by 
simulations and theory. We set the interaction between B patches to zero and calculate the 
phase diagram, as the ratio of the AB and the AA interactions, e^^, varies. We find that, in 
the same range of parameters as in 3D off-lattice models, the liquid-vapor diagram exhibits 
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a re-entrant or "pinched" shape, and there is an empty fluid regime. In addition, below 
^*AB — \ condensation ceases to exist, and the re-entrant regime disappears for > |, 



m 
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ine with the results for off-lattice 3D models and the predictions of Wertheim's theory 



13| . This suggests that the thresholds predicted by Wertheim's theory are exact and 



universal, i.e. independent of the dimensionality of the system and of the lattice structure. 
Finally, the 2A2B model exhibits an order-disorder (0-D) transition that is ferromagnetic 
for l<e\B< i 

The use of 2D lattice models allows the simulation of larger systems enabling us to 
establish the nature of the critical points and to investigate the structure of the network 
liquid phase in the empty fluid regime, where the size of the "voids" increases rapidly as 
the temperature decreases. We find that the liquid-vapor critical points are in the 2D Ising 
universality class, with a scaling region that decreases as the temperature decreases. The 
simulation results also indicate that the line of 0-D transitions intersects the condensation 
line at zero temperature and density, at a multicritical point, or at a very low temperature, 
at a critical end-point. The analysis of this region requires the simulation of larger systems at 
extremely low temperatures, which becomes prohibitive even for 2D patchy particle lattice 
models. 

In order to proceed we consider a low-temperature-model (LTM) that describes the 
asymptotic behavior of 2D patchy particle models at low temperatures, and use a cluster 
algorithm that enables the efficient simulation of these low temperature systems. Finally, 
we derive asymptotic expressions based on Wertheim's theory, for the liquid branch of the 
binodal and the 0-D transition that suggest, in line with the simulation results, that the 
transition lines meet, at a multicritical point, at zero temperature and density. 

The paper is arranged as follows: In section II we describe the patchy particle model, the 
mapping of the full lattice limit and the simulation methods. In section III we present the 
results for the global phase diagram of a system with a re-entrant binodal. We compute the 
binodals, analyze the nature of the liquid-vapor critical points, and discuss the topology of 
the global phase diagram for systems with a ferromagnetic ordering transition (re-entrant 
regime). In section IV we introduce the LTM and the simulation techniques developed to 
sample low temperatures efficiently. We compute the binodals and the ferromagnetic or- 
dering transition and discuss the topology of the global phase diagram. Then, in section V 
we address the zero temperature and zero density limit theoretically. We derive asymptotic 
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expressions for the condensation and 0-D transitions based on Wertheim's theory for asso- 
ciating hquids, in the hmit of strong AA bonding. We conclude, based on the asymptotic 
analysis, that the condensation and 0-D lines meet at a multicritical point, at zero temper- 
ature and density. In section VI we make some concluding remarks and in the Appendix 
provide details of the calculation of the starting point of the liquid-vapor equilibrium of the 
LTM, used in the Gibbs-Duhem integration of the liquid branch of the binodals. 



II. THE 2A2B MODEL 



The model consists of particles with four patches, two of type A and two of type B, on 
a square lattice. The lattice sites are either empty or occupied by one single particle. The 
patches A and B are aligned along one of the two lattice directions (See Figured]). There 
are two configurations for each occupied site: (1) A patches aligned along ±x and B patches 
aligned along ±y , and the symmetric configuration with (2) A patches aligned along ±y 
and B patches aligned along ±x. The potential energy, U, is the sum of pair interactions 
between nearest-neighbor (NN) particles on the lattice and is written as: 

W = -eAA-f^AA - (^AB-f^AB " ^BB-f^BB] (1) 

where A/'a^ is the number of a(3 bonds, i.e. lattice bonds between NN occupied sites con- 
necting patches a and (3. 

This model is a lattice realization of the patchy particle models with distinct patches 



introduced in 



5| and investigated in the context of empty network fluids [ll|, [12] • In 
line with previous work we take caa = e as the energy scale (e > 0), and focus on systems 
where the B patches do not interact, i.e. ebb = 0. The interaction between A and B 
patches varies although most of the results are for systems in the re-entrant regime, i.e. 
< tAB < e/2. Taking into account that each particle carries two A patches and that a 
patch can participate, at most, in one bond we can write: 

2N = 2AfAA+AfAB+AfAo; (2) 

where N is the number of particles in the system and A/^o is the number of patches that are 
not bonded. Combining Eqs. ([1]) and ([2]) for e^s = we get: 

U/e = -N + Mab - ^) + ^-^^0. (3) 



From Eq. (|3]) It follows that AA bonds are favored when cab < e/2, while AB bonds 
are favored when e^s > e/2. At low temperature, most of the A patches are bonded and 
the network fluid consists of AA chains connected by a small number of AB branches in 
the former case while the network is almost fully branched in the latter. The special case 
^AB = corresponds to a self-assembling rigid rod (SARR) model that was studied on 2D 



lattices recently 



21 



2q . In the SARR model, a continuous 0-D transition is found to be 



the only feature of the phase diagram. At low temperatures the particles form long rigid 
rods, through AA bonds, which undergo an orientational ordering transition, in the 2D Ising 



class on the square lattice and in the q=3 Potts class on the triangular one 



'22 



21. The 



SARR model has no liquid- vapor transition as adjacent rods do not interact energetically. 



The full lattice limit 



The full lattice limit of the SARR model on the square lattice has been mapped on 



to the Ising model 22, l26|. This was achieved by establishing a correspondence between 
the particle orientations of the SARR model and the spins ±1 of the Ising model. The 
total energy of both models is then computed by adding the contributions of elementary 
plaquettes, consisting of a square with four sites enclosing an elementary lattice cell. The 
mapping between the Ising and the full lattice limit of the SARR model, is established for 
any plaquette configuration and the critical temperature of the model is identified with the 
exact result of the critical temperature of the corresponding Ising model j^]- Following this 
procedure a similar mapping is established for the 2A2B model. At full lattice occupancy, 
the 2A2B model undergoes an Ising 0-D transition, at the reduced temperature 22|: 



(4) 



2eln {1 + V2) 

where ks is Boltzmann's constant. The ordered phase is stable at T < Tc. When esB = 0, 
the ordered phase is ferromagnetic (particles aligned in the same direction) if e^B < e/2, 
and antiferromagnetic otherwise, eA_B > e/2. 

Note that when eBB = and bab = e/2 there is no 0-D transition. Inspection of 
Eq. ([3]) reveals that at full lattice occupancy every A patch is bonded, and therefore all 
the configurations have the same potential energy. This degeneracy does not hold when 
vacancies (empty sites) are present but the free energy is still dominated by entropic terms 




FIG. 1: Illustration of the model. Top: One particle with 2 A and 2B bonding sites or patches with 
the A patches aligned along x and the B patches aligned along y. Bottom: Two particles forming 
a A A bond along x. 

that prevent the system from ordering. 



B. Simulation methods 

We aim at computing the global phase diag ram o 
Carlo Simulation. Based on previous results 
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12 



22 



the 2A2B model through Monte 



23| a low temperature critical line 



corresponding to the 0-D transition, the locus of which at p = 1 is known exactly through the 
mapping to the Ising model is expected to occur; in addition, a liquid- vapor first-order 



8 



transition ending at a critical point (for cab above a certain threshold) is also expected. 



The 0-D transition is located using techniques analogous to those described in |23|. We 
fix the temperature and system size and, by means of the simulated tempering algorithm, 
compute the properties for different values of the chemical potential around the expected 
critical point. By using appropriate finite-size scaling analysis we obtain estimates of the 
critical parameters, fic{T) and Pc{T), in the thermodynamic limit. 

The liquid vapor equilibria (LVE) is computed using a combination of Wang-Landau mul- 
ticanonical simulation (WLMC) [27| and Gibbs-Duhem integration (GDI) procedures [28 1 



adapted to lattice models [29|, 



30| . The WLMC methodology was described previously, in- 



cluding the details specific to lattice models 



29 



31 



32| . WLMC simulations, combined with 



finite-size analysis techniques, are very efficient in locating the liquid- vapor critical point, 
and in computing the phase diagram at temperatures not far from it. At low temperatures 
we found it useful to resort to GDI schemes. 

We run WLMC simulations and locate the LVE at a given temperature and system size 
by searching for the value of the chemical potential, /io(L,T), that maximizes the density 
fiuctuations: Sp = [< > — < p >^]^^^. Under these conditions we compute the average 
density, pm = Pm{L, T, /iq), and the moments of the density distribution, nik =< {p—pm)^ >, 
in order to calculate the ratio = m/^/m^, which is related to the fourth-order Binder 



cumulant 



33l | . We establish whether at the chosen temperature T there is LVE by analyzing 



the dependence of on the system size. LVE occurs below the critical temperature, where 
at po the density distribution function exhibits two peaks that become sharper as the system 
size L increases. This implies that gi{L) decreases as L increases and approaches g^ = 1 
in the thermodynamic limit. Above T^, gi{L) increases with L and approaches 5'4 = 3 
(Gaussian distribution) in the thermodynamic limit. At the critical temperature, finite-size 



scaling arguments [27, 



34 



35], indicate that (for sufficiently large systems) g^^L) takes a 



non-trivial value that depends on the boundary conditions and on the universality class of 
the transition. 

We estimate system size dependent pseudo-critical points: \Tc{L) , pc{L)] by imposing 
that (y'4(L, T) takes the value corresponding to the 2D I sing universality class [36l|. Numerical 



details of these calculations may be found elsewhere [27|, l32|. We proceed to estimate the 
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critical temperature and density in thermodynamic limit, using the scaling equations: js^ 

p,(L) - oc (5) 

T,(L) - T, oc L-'^-\ (6) 

where v is the correlation length critical exponent {u = 1 for the 2D Ising class); and 
A = 6/ V, where 9 is the correction to scaling critical exponent. There is some controversy 



37 



41| concerning the value of A for systems in the 2D Ising class, as a number of simple 



models (e.g. 2D Ising) have no irrelevant operators [36|, |37[. One then expects, A = 4/3 



4l| in general or A = 7/4 36|, |38|] in the absence of irrelevant operators. Taking this into 



account, we computed three estimates of the critical temperature (the same scheme applies 
to the critical chemical potential), using: A = 4/3; A = 7/4; and considering A as a fitting 
parameter. 

At temperatures below we fit the system-size dependent LVE results to the scaling 
equations: 

x(L,T) -x(T) oc L"^; (7) 

where x(L, T) is the finite-size result for the property x, and x{T) (obtained from the fit) is 
the estimate of the property in the thermodynamic limit; d = 2 is the spatial dimensionality 
of the system. We have obtained very precise values of the chemical potential at coexis- 
tence, which were subsequently used in the GDI to compute the LVE in a wider range of 
temperatures (away from the critical point). Within the GDI we run sequences of two phase 
(liquid and vapor) simulations using larger system sizes (than those feasible with WLMC) 
allowing us to sample lower temperatures, in the empty liquid regime. 

In the computation of the critical parameters, described above, we assumed that the 
critical point of the LVE is in the 2D Ising universality class. As this is not yet established, we 
proceed to analyze the scaling behavior of the pseudocritical parameters, and the moments of 
the density distribution, Pl{p), at the pseudocritical points. The finite-size scaling behavior 
of Spc{L) satisfies: [42 ] 

L^'/^SpciL) ^ao + aiL-^; (8) 

where /3' is the critical exponent of the order parameter (/?' = 1/8 for 2D Ising). We 
expect the shape of the critical density distribution, Pl{p), to approach that of the crit- 
ical Ising 2D magnetization P^^*"^(A^), for large system sizes [34.]. Deviations occur for 
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small systems due to corrections to scaling associated to irrelevant fields and field-mixing 



contributions 



34| . Thus, in addition to checking the scaling of Tc(L), Pc{L) and 6pc{L), 



we compare the asymptotic values of the reduced moments of the density distribution, 
(75 = m^jmcl'^ and = m^/m^^ to the critical Ising 2D values, (75 = and ^ 1.4556 



III. RESULTS FOR THE 2A2B MODEL 



We start by illustrating, in Figure [21 typical configurations of the coexisting phases for a 
system in the re-entrant regime, with e^^ = eAs/^AA = 0.40, at three temperatures. We note 
that the density of the liquid decreases rapidly as the temperature decreases. The (network) 
liquid phase is characterized by voids (regions without particles) that increase as the tem- 
perature decreases. This observation implies that larger system sizes are required at lower 
temperatures, in order to sample adequately the increasing length scales that characterize 
the empty hquid phase. 



A. The phase diagram 

We consider 2A2B models characterized by different e^^. After preliminary WLMC tests 
we choose appropriate subcritical temperatures and compute the liquid-vapor equilibria by 
extrapolating to the thermodynamic limit the results of several system sizes. We then select 
a temperature (for each model) as the starting point of the GDI. These are collected in 
Table [T] 

In Figure [3] we plot simulation and theoretical results for the liquid-vapor binodal of 
the 2A2B model with e^^ = 0.40. The binodal has the "pinched" or re-entrant shape, 
characteristic of 3D off-lattice patchy particle models, with two A patches and | < e^^ < ^ 
jnl . \l2 \. The coexisting liquid density vanishes rapidly as the temperature decreases and 
the model exhibits an empty fluid regime. The theory (based on Wertheim's theory for 
associating fluids discussed in section V) describes the re-entrant behavior of the binodal 
and gives a good estimate of the critical temperature but underestimates the coexisting liquid 



density, as in related 3D off-lattice models 



11 



15| . The cornputed percolation threshold, 



for clusters of particles connected by bonds between patches [l2|], is also shown in Figure |3l 



The simulation results suggest that the percolation line intersects the LVE binodal at the 
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(a) e*ab=0.4, r=0. 1 2, Vapor Phase (b) eab=0.4, T=0. 1 2, Liquid Phase 
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(c) e*ab=0.4, T*=0.08, Liquid Phase 
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(d) e*AB=0.4, T*=0.06, Liquid Phase 
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FIG. 2: Representative configurations of the 2A2B model with = 0.40 and L = 128 at liquid- 
vapor coexistence, and several reduced temperatures T* = ksT/e. (a) Vapor phase at T* = 0.12; 
(b) Liquid phase at T* = 0.12; (c) Liquid phase at T* = 0.08; and (d) Liquid phase at T* = 0.06. 
Particles are represented as segments of unit length oriented in the direction of the A patches. 
Note that the liquid becomes emptier as the temperature decreases. 



critical point, in line with results for the 2D Ising model 43|. This contrasts with the results 
of Wertheim's theory (details of the theoretical methodology can be found in Ref. Il3l ) and 
the simulation results of 3D off-lattice models, where the percolation line intersects the LVE 
binodal on the vapor side jisj ]. 



B. The nature of the liquid- vapor critical points 

In Figure m we illustrate the scaling behavior of the critical parameters and the moments 
of the density distribution function, at the LVE pseudo-critical points, with the system 
size, for two 2A2B models with: = 0.40 and = 0.50. The observed behavior is 
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^AB 




L(GDI-LT) 


L(GDI-HT) 


0.375 0.10 


-1.01582(2) 0.155(1) 0.146(3) 


512 


512 


0.400 0.13 


-1.04302(2) 0.249(1) 0.226(1) 


512 


512 


0.450 0.15 


-1.08487(2) 0.3520(2) 0.3376(2) 


512 


256 


0.500 0.20 


-1.16526(2) 0.3891(4) 0.3242(5) 


128 


64 


0.550 0.22 


-1.22628(2) 0.4251(2) 0.37544(3) 


128 


128 


0.600 0.25 


-1.30164(2) 0.4426(3) 0.3704(2) 


128 


128 



TABLE I: Liquid vapor equilibria from WLMC simulations. The temperatures are those used as the 
starting points of subsequent GDIs. L(GDI-LT) and L(GDI-HT) correspond to the largest system 
sizes used in the GDI for temperatures below and above (respectively) the starting temperature Tq. 
Error bars, between parentheses, are given in units of the last digit and correspond to a confidence 
level of about 95 %. 




panel: Results of Wertheim's theory. See the legends for details. 

consistent with criticality in the 2D Ising universality class. Note, however, that the system 
size dependence of the critical properties is stronger in the system with e^^ = 0.40. 

The results for the critical temperature and the critical chemical potential hardly depend 
on whether we use A = 4/3; A = 7/4 or A being considered as a fitting parameter (See Eql6]). 
In the latest case the effective values of A are always larger than A = 1 (See the effective 
values Xeff in Table HTll . In particular, for the largest values of the effective values of A 
are consistent with A = 4/3. For = 0.35,0.40 the uncertainty in the effective value of 
A is too large to discriminate between the two scaling scenarios. Nevertheless, the fact that 
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0.02 0.04 0.02 0.04 ' 0.02 0.04 

1/L 1/L 1/L 



FIG. 4: Scaling of the critical parameters and of the moments of the density distribution function, 
at the pseudo-critical point, with the system size, for two 2A2B models; one model is in the re- 
entrant regime and the other is at the boundary to normal liquid behavior. The scaling results 
are fully consistent with 2D Ising criticality. The pseudocritical parameters, Pc{L), Tc{L) and 
5pc{L), follow the 2D Ising scaling laws, Eqs. ([5]), ^ and ([8]); while the moment ratios, g^, and 
qq, approach the 2D Ising values as L increases. It is also clear that the scaling region decreases 
as decreases. 

the effective values of A satisfy A > 1, supports the hypothesis that the LVE critical point 
of the 2A2B model is in the 2D universality class. 

In Figure O we plot the liquid- vapor phase diagram for various 2A2B models in the re- 
entrant and normal regimes. Numerical results for the critical points are collected in Table 
im The results shown for and fic are those extracted from the fitting scheme with fixed A 
that provides the best agreement with simulation data (A = 7/4 for < 0.40, and A = 4/3 
for > 0.40). As expected both the critical temperature and the critical density increase 
with e\^. A significant change in the binodal, however, occurs at = 0.5. For models 
with e^^ > 0.5 the liquid density approaches p = 1 as the temperature vanishes; while for 
models with e^^ < 0.5 the liquid density decreases with temperature, at low temperatures 
and seems to approach p = as the temperature vanishes. This conclusion is based on 
theoretical results (see section V) and confirmed by computer simulations of 3D off-lattice 
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FIG. 5: Simulation results for the liquid-vapor binodals of different 2A2B models in the re-entrant 
and normal regimes. Prom top to bottom: = 0.60, 0.55, 0.50, 0.45, 0.40, 0.375. 



models 
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12| . The simulation of systems at vanishingly low temperatures is hindered 
by two factors: the usual problems of sampling at low temperatures, and the emergence of 
diverging length scales, namely the size of the voids in the empty liquid phase. 

The model with e^^ = 0.50 exhibits an intermediate behavior. Simulation results suggest 
that at T = the density of the liquid phase at equilibrium with the vapor approaches a 
finite density, p ~ 0.87. Note also the dashed lines that continue the hquid branches of the 
models with e^^ < 0.50. These lines were computed using a related model that captures the 
phase behavior of the 2A2B patchy particle models at low temperatures, to be described in 
section IV. 

A final important question remains. The 2A2B patchy particle models undergo, in gen- 
eral, two thermodynamic transitions, a first-order liquid-vapor transition at low densities; 
and a continuous 0-D transition at high densities (See the phase diagram for e^^ = 0.40 in 
Figure |3]) . Previous results [Ul, Il2| and those discussed here suggest that for appropriate 



values of e^^ both the vapor and the liquid branches of the LVE approach zero density at 
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23| also indicate 



zero temperature; on the other hand the results for the SARR model 
that the 0-D transition approaches zero density at zero temperature. 

Several scenarios are then plausible for the global phase diagram of the 2A2B patchy 
particle model depending on where the 0-D critical line intersects the LVE line; this can 
happen at a finite temperature (either at a critical end point or at a tricritical point) or at 
T 



0. The simulation results for e^^ 



0.45, and e^^ = 0.40, discard the possibihty of an 
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AB 


C 


pt Pc Kff 




-^max 


0.375 


0.1160(10) 


-1.0251(6) 0.170(5) 2.3± 1.6 


64 


128 


0.400 


0.1402(2) 


-1.0509(2) 0.242(1) 1.7± 0.3 


36 


160 


0.450 


0.1768(2) 


-1.1087(2) 0.3279(5) 1.3 ± 0.3 


20 


64 


0.500 


0.2068(1) 


-1.1715(1) 0.3782(3) 1.3 ± 0.2 


12 


72 


0.550 


0.2340(1) 


-1.2387(1) 0.4109(4) 1.2 ± 0.3 


10 


56 


0.600 


0.2599(2) 


-1.3099(2) 0.4339(3) 1.3 ± 0.2 


8 
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TABLE II: Numerical results for the liquid-vapor critical points of different 2A2B models in the 
re-entrant and normal regimes. The effective exponents Ae// are obtained from the fits of the 
pseudocritical temperatures using Eq. ([6]) with v = 1. Error bars correspond to a confidence level 
of about 95%. In the cases of Tc and the error bars extend over the results of the three fitting 
schemes (A = 4/3, A = 7/4; and A as adjustable parameter). 



upper tricritical point, as the critical temperature of the LVE is higher than the exact result 



for the 0-D transition at p = 1 and previous 



22| and present results (See Figure |3]) suggest 



that the critical temperature of the 0-D transition increases with the density. 

The question remains whether the transitions meet at T = or at a finite temperature 
critical end point. For the models with e^^ < 1/2 considered here we computed the order 
parameter on the liquid branch as obtained from the GDI, for several system sizes, and 
found that the network liquid phase at LVE is orientationally disordered. If a critical end 
point exists then it has to occur at lower temperatures than those accessible by simulations 
of the 2A2B patchy particle model. To proceed we consider a related model in the next 
section. 



IV. PHASE BEHAVIOR AT LOW TEMPERATURES 

Let us consider 2A2B models with ess = 0. Defining A = | — we re- write Eq. ([3]) 

as: 

U/e = -N + XArAB + lAfAo, (9) 

with A = for e^^ = e/2, and A = 1/2 for = 0. Now consider < A << 1/2. 
At sufficiently low temperatures A/^o (the number of non-bonded A patches) is negligible 
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with respect to both Mab and XMab, and the thermodynamics of the model is determined 
by A* = Xe/ksT. This suggests that the LVE of 2A2B patchy particle models at low 
temperatures, and small A, may be collapsed (approximately) onto a single curve. 



A. The Low Temperature Model 

The previous discussion suggests the consideration of a related low temperature model 
(LTM) with interaction energies: uaq = oo (i.e. non-bonded A patches are disallowed), 
uaa = ubb = ubo = 0, and uab = Ae. 

In order to compute the LVE we assume that the vapor phase at low temperatures has 
zero density. This results from the fact that, as all A patches are bonded, the particles 
must belong to a network that percolates in, at least, one direction; since the vapor does 
not percolate its density must vanish in the thermodynamic limit. We compute the liquid 



branch of the LTM using Gibbs-Duhem integration, with t 



lermodynamic variables A* and 



[Pfi) = fi/ksT. The differential equation to be solved is:|30| 

Ndi^jj) - AfABdX* = (10) 

V dX* <N> ^ > 

where N and Mab correspond to the liquid branch (the vapor has zero density). The starting 

point for the integration is (/3/i) = (/3/i)o, A* = 0. The calculation of (/3/i)o is discussed in the 

Appendix. At full lattice occupancy the LTM is equivalent to the original patchy particle 

model, since all A patches are bonded, and thus the LTM may also be used to compute the 

0-D transition. 

In the LTM every A patch is bonded (either to another A or to a 5 patch). In addition, 
the density of the liquid phase decreases rapidly as A* increases, in the empty fluid regime. 
A standard algorithm involving single particle moves is useless under these conditions. In 
order to sample the LTM we develop an efficient cluster algorithm that is described below. 



B. The cluster algorithm 

The LTM algorithm is based on three types of moves: (a) Rotation of particles (only if 
the four NNs of the particle are occupied); (b) Insertion of a sequence of aligned particles; 
(c) Deletion of a sequence of aligned particles. 
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The rotation move is straightforward. One particle with four NNs occupied is selected at 
random (if there is any), and then one of its two orientations is chosen with a probability 
proportional to its Boltzmann factor. 

The insertion / deletion of a sequence of aligned particles is carried out as follows: A lat- 
tice site is chosen at random: If the site is occupied then a deletion attempt is performed. It 
starts by identifying the linear cluster of particles linked to the selected one by an unbroken 
sequence of AA bonds. Such cluster either percolates through the periodic boundary condi- 
tions (PBC) or ends at two AB bonds. If the removal of the cluster leads to an unbonded 
A patch the deletion attempt is rejected, otherwise the acceptance criterion (defined below) 
is applied. 

If the chosen site is empty then one direction, s = 1,2, is chosen at random. A linear 
AA cluster of occupied sites is built along the chosen direction (on both sides), the bonding 
criterion being that the NN position is empty. The process stops when the cluster percolates 
through the PBC or when occupied sites are found at both ends. The acceptance criterion 
(defined below) is then apphed. 

It is straightforward to compute the change in energy when inserting or deleting a LTM 
cluster. The cluster either percolates through a sequence of AA bonds {AU* = 0) or 
terminates at both ends with AB bonds: AU = ±2A. Considering that positions (not 
insertions/deletions) are selected at random, the acceptance probabilities are: 

A{N\N + AN) H \ k j 

where AU = U^^/^jy — U^, and the factor 2 arises from the two orientations of the inserted 
cluster, of length AA^ lattice sites. 

C. Simulation Results 

The GDI requires as input a point on the LVE binodal, which was taken to be (/3/i)o, the 
reduced chemical potential at A* = 0. The chemical potential at zero pressure (the vapor 
phase has zero density) is obtained via thermodynamic integration {44 1 from {13 fx — > 00), as 
the partition function in the full lattice limit and A = is known exactly: Q = 2^ . We 
carried out the calculation for different system sizes L = 16, 32, 64, ■ ■ ■ , 256 and found that 
the size dependence of /i is negligible, obtaining {f3fi)o = —0.78940(2). This is consistent 
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0.2 0.4 0.6 0.8 



P 

FIG. 6: Phase diagram of the LTM (LVE and order-disorder transition); re-scaled hquid- vapor 
binodals of the 2A2B patchy particle models with < 0.50, and order-disorder transition for 
^*AB ~ 0-^0 (Symbols are explained in the legends). The liquid branch of the LTM model is 
computed for a system with L = 2048. 

with the estimate from the GDI of the patchy particle model with = 0.50, which 
gives /i/e — 1 — 0.789T*, at low temperatures. The LVE is obtained using Gibbs-Duhem 
integration. Several system sizes are considered to check the system size dependence of 
the results. In Figure E] we test the accuracy of the LTM to describe the coexisting liquid 
densities of the 2A2B patchy particle models, at low temperatures. Clearly, the LTM results 
converge to those of the patchy particle models as the temperature is lowered. We found 
that as A* increases (scaled temperature decreases) larger systems are required to obtain 
consistent results for different system sizes (as was observed in the simulation of the 2A2B 
patchy particle models). Indeed, the line corresponding to the LTM liquid branch in Figure 
E] is plotted for scaled temperatures higher than those where the results for systems with 
L = 1024 and L = 2048 start to show significant differences: i.e. A* ~ 1.88. This is due to 
the rapid increase of the size of the voids in the empty liquid at low temperatures, which 
hinders the simulations of the LTM at larger A* due to the system size requirements and 
the loss of efficiency of the simulation algorithm. 

Finally, we have computed the 0-D transition of the LTM, which is almost indistinguish- 
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able from that of the 2A2B model with e^^ = 0.40, after proper re-scaling (both are plotted 
in Fig. |6]). 

Despite the difficulties in simulating the LVE of the LTM when A* > 1.9, the numerical 
results suggest that the ratio pl{X*) / pooi^*) decreases with A*, for A* > 1.60. Considering 
that the LTM describes accurately the low temperature phase diagram of the 2A2B patchy 
particle models, we conclude that the most likely topology of the phase diagram of this 
class of models is characterized by a multicritical point at T = and p = 0, where the 
liquid-vapor and the 0-D transitions merge. 



V. TOPOLOGY OF THE PHASE DIAGRAM: THEORY 



A. Wertheim's Theory 

In this section we address the topology of the phase diagram by resorting to theoreti- 
cal/analytical techniques. The thermodynamics of the 2A2B patchy particle model can be 
described using Wertheim's first order perturbation theory (WPT), which accounts accu- 
rately, in the low density limit, for the effect of association 11, 12|. The reference free energy 
Fref is that of an ideal lattice gas, 

^ = lnp+l^ln(l-p), (13) 

where p is the density, 
given, within WPT, by 



"?he perturbation term Fj, includes the bonding contribution and is 



^ = 2lnXA-XA + 2lnXB- Xb, (14) 

where X^ is the probability that a bonding site of type a is unbonded. These probabilities are 
related to the thermodynamic quantities through the laws of mass action (i.e. by considering 
bond formation as an equilibrium chemical reaction), which are, for particles with 2 A and 
2B bonding sites, 

X, + 2p A,„X2 + 2p A,^X„X^ = 1 , (15) 
with a = A,B and P a. The quantities A^^ are given by, 

= Vai3 [exp(/3e„/3) - 1] , (16) 
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0.00, 




FIG. 7: Phase dia gram of the 2A2B model, based on Wertheim's theory, for several values of e^^- 



with Vai3 the volume of the a/3 bond (in units of the volume of a lattice site). Note that, 
since ess = 0, A^^ = 0. We take Vb = 1/2, in order to maximize Vb while disallowing more 
than one bond between two patches and more than two patches per bond (see Figured]). 

Using equations ( fT5|) in ( IT4l) . Fb is obtained as a function of p and T. From the Helmholtz 
free energy F = F^ef + Fb-, one can obtain the pressure and the chemical potential and 
calculate the phase diagram. Figure [7] shows the results of this calculation for several values 
of e\Q. Comparison with the results of simulations (section III) reveals that the theory 
describes correctly the re-entrance of the liquid branch for e^B < 0.5e and the constant 
density of the liquid branch at low temperatures when eAB = 0.5e. The theory also predicts 



{4, luj that no liquid- vapor coexistence occurs when eAB < e/3, in line with the results of 
ihe simulations for these values of the parameters. Note that, as in 3D off-lattice models 



11 



, Il3| , there is almost quantitative agreement between the critical temperatures obtained 
by theory and simulations, while the theory underestimates systematically the density of 
the coexisting liquid branch. 



B. The liquid branch of the binodal and the order-disorder transition 



Wertheim's theory, as described in the previous section, is not capable of describing the 
0-D transition. The ordering is driven by the excluded volume of the chains formed at low 
temperatures 21|, and this effect is not included in ( IT3ll nor in (IHj). Based on previous 



works 



11 



21| . we proceed to derive asymptotic expressions for the 0-D transition and the 



liquid branch of the binodal. This analysis shows that the line of LVE is not intersected by 
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asymptotic 







□ 




o 


V=AA="-375 



0.25 



FIG. 8: Binodals for < 0.5 at low densities and temperatures calculated using Wertheim's 
theory. The black line is the asymptotic result (jlSp for the liquid branch of the binodal. Note that 
the temperature of each binodal is rescaled by A = 1/2 — e^^. 

the 0-D line at any finite temperature, and thus a critical end point does not occur in this 
model. 

The asymptotic limit for the liquid branch of the binodal of 2A2B models in the re- 



entrant regime, e^s < 0.5e, is obtained using the results of Taking into account that 
the reference free energy is given by ( !T3|) . the asymptotic pressure, i.e. the pressure in the 
limit of strong AA association within WPT, is given by, 

^2 



/3p = aop2 - aip2 + 



El 
2 ' 



(17) 



with tto = {2Aaa)~^ and ai = 2AabO'0- The first term vanishes when all the A patches are 
bonded (see section IV) and under these conditions p ~ at coexistence. The coexisting 
liquid density, pe, is then approximated, at low densities and temperatures, by, 

^ = ^ . (18) 

eAA^ lnp^-21n2 ^ ^ 

The asymptotic liquid density is plotted in Figure [8] together with the binodals of various 
2A2B models, obtained using Wertheim's theory. It is clear that Eq.f llSI) is the asymptotic 
limit of the liquid branch of 2A2B patchy particle binodals, at low densities and tempera- 
tures. 

""inally, we turn our attention to the 0-D transition of the 2A2B patchy particle models. 



In 



2l[ | the SARR model, which is the limit of the 2A2B model when eAB = (i.e. when 



only A A chains are formed) , was investigated and the contribution of the excluded volume 
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FIG. 9: Asymptotic results for the liquid branch of the binodal (full line) given by (jl8l) and for the 
order disorder transition (dashed line) given by (j26j) . 

of two chains was included in the free energy via an Onsager like approximation. An order 
parameter A = Px—Py (with pa being the number density of particles with A patches aligned 
along a) was defined, and it was shown that the field ho associated to A, is, 



In 



A 



with X and Y given by. 



Px = exp(-/3e) 
Py = exp(-/3e) 



X 



;i -x)2^ 

Y 



(19) 

(20) 
(21) 



(1-F)2- 

The 0-D line is found by solving (^)^^o = 0. 

The 2A2B model differs from the SARR model by allowing the formation of AB bonds. 
In the limit of strong AA bonds the field h, conjugated to the order parameter A, may be 
taken to be. 



I3h = l3ho + 



AB 



OA 



(22) 



where Jab is the contribution of the AB bonds to the free energy density. If this free energy 



is calculated within WPT and taken in the limit of strong AA bonds 

l3fAB = -4^201 (^plpy + pIpx) , 
one obtains a simple expression for h, 
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(3h = (3ho - 2V2ai Q(px^ 



Py + Py ^Px) + Px + Py 



(23) 



(24) 
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The 0-D line is now calculated by solving (^)^_q = 0, which, in the limit of low densities 
and temperatures, yields, 

ooP^ - + 5aipi = 0. (25) 

In line with the derivation of f lTS]) . we neglect the first term of Eq. fl2^ . and obtain for the 
0-D transition line, 

c^aA lnp-51n2' ^ ' 

The asymptotic liquid binodal ( ITSj) and the asymptotic 0-D line ( 126|) are plotted in Figure 

[9l they do not intersect at finite temperature and thus the global phase diagram of the 

2A2B model does not have a critical end point. In other words, the empty liquid regime is 

not pre-empted by the ordered (liquid) phase. 



VI. CONCLUSIONS 



We investigated a simple patchy particle lattice model consisting of particles with four 
bonding sites, two of type A and two of type B, on the square lattice, and computed its 
global phase behavior by simulations and theory. We have set the interaction between B 
patches to zero and calculated the phase diagram as the ratio between the AB and the AA 
interactions, e^^, varies. In line with previous work on three-dimensional off-lattice models, 
we have shown that the liquid- vapor phase diagram exhibits a re-entrant or "pinched" shape, 
for an identical range of e^^, suggesting somewhat surprisingly that this range - and the 
corresponding empty fluid regime - is independent of the dimensionality of the system and 
of the lattice structure. 

In addition, the use of low- dimensional lattice systems allowed the simulations of much 
larger systems enabling us to establish the nature of the liquid-vapor critical points, which 
were found to be in the Ising 2D class, both in the re-entrant and the normal liquid models. 
While in the normal liquid regime the scaling regions are typical of models in the 2D Ising 
universality class, in the re-entrant liquid regime the scaling region decreases rapidly as the 
critical temperature (or e^^) decreases. Our theoretical and simulation results also suggest 
that the Ising scaling region vanishes as the critical temperature vanishes, in line with the 
presence of a multicritical point at zero density and temperature. 

The patchy particle models on the square lattice exhibit an 0-D transition at fixed 
density, as the temperature is lowered below the bonding temperature. This transition 
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is anti-ferromagnetic for normal liquid models and ferromagnetic for models with a re- 
entrant liquid regime. In the latter models, the results of simulations of an appropriate 
low-temperature-model that describes the asymptotics of the particle patchy systems at low 
temperatures, together with an efficient sampling cluster algorithm, suggest that the line 
of 0-D transitions intersects the condensation line at zero temperature and zero density. 
This topology of the phase diagram is corroborated by an asymptotic theoretical analysis of 
the hquid branch of the binodal and of the 0-D transition, based on Wertheim's theory for 
associating fluids. The theory is exact at zero density, lending support to the results of the 
asymptotic analysis in the low temperature, low density region. 

In summary, the results of simulations and of theoretical analysis strongly suggest that 
the line of 0-D transitions intersects the condensation line at a multicritical point at zero 
temperature and density, for patchy particle lattice models in the re-entrant liquid regime. 
The global phase diagram of off-lattice patchy particle models, in 2D and 3D, is further 
complicated by the presence of stable solid phases. These phases may pre-empt the empty 
fluid regime rendering the zero temperature zero density multicritical point metastable. 
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VII. APPENDIX: COMPUTING THE INITIAL POINT FOR THE GDI OF THE 
LTM 

The hmit A* = of the LTM is an athermal model, where all allowed configurations have 
zero potential energy. A configuration is allowed if (and only if) every patch of type A is 
bonded. In this model a first-order transition, corresponding to the transition of the 2A2B 
model with cab = e/2 at T = 0, is expected to occur. The results of a series of simulations 
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of the LTM with increasing / decreasing values oi (3/1 indicate that the transition at A = 
is indeed first order. The coexisting vapor phase is found to have vanishingly small density, 
in line with the results for the 2A2B model. 

The value of at the transition is computed using thermodynamic integration, as the 
partition function for A* = in the full lattice limit is known exactly: 

Q{N = M,M)=q^, (27) 

where q is the number of particle orientations, q' = 2 for the square lattice. If the number of 

vacancies is small: M — N « M, we can write an approximate expression for the partition 
function, by assuming that the number of NN pairs of vacancies is negligible: 

QIN, M) ~ QJN, M) = 2"-''<"-"> ( ] = 2™-"" | | . (28) 

The factor 2~^(^~-^) arises as an isolated vacancy ehminates the possibihty of having two 
different states (orientations) at the vacant site and at its four NNs (that are assumed to be 
occupied), since in the LTM orientations with A patches pointing to the vacant site are not 
allowed. 

We can assume the vapor phase density to vanish, and thus the pressure at coexistence 
also vanishes. As the partition function at full coverage is known, we can compute the 
equation of state of the high density phase and the value of the chemical potential at the 
transition ni^, which satisfies: /3p{/3fj,iy) — 0. 

In order to derive a procedure to compute /3/ii„ we consider, first, the approximate par- 
tition function Qo{N, M). The corresponding grand canonical partition function is: 

go(/3//,M)) = 2-^^ J] exp[(/3// + 51n2)iV]; (29) 



N=0 

which can be summed to give: 



N 



go(/3/x, M) = 2-^^ [1 + exp ■ (30) 

where j3^' — (3/1 + 51n2. The pressure, at this level of approximation, is written: 

(/?/.', M) = ^lnQo(/3/^,M) = -41n2 + In (l + e'^'^') ; (31) 
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At the same level of approximation, the density is easily computed: 



E 



M 

N=l 



M 
N 



M 



E 



M 

N=l 



M 
N 



1 + e'^^'' 



(32) 



We now define the fugacity fraction (f) as: 



(33) 



^ 1 + e^^'' ' 
we obtain, 

pW(0) = 0, (34) 

(0) = -4 In 2 - In (1 - 0) . (35) 

In the Grand Canonical Ensemble for processes at constant temperature and constant vol- 
ume, we have: 

d{(3p) = pd{(3fi), (36) 
which can be integrated to give the pressure as a function of the fugacity fraction: 



/3p(0) = /3p(°)(0)+ / d0i[p(0i)-0i 



pp{(j)) = -41n2 - In (1 - 0) + j uc/^i 



p(0i) - 01 



(37) 



(38) 



The integrand in Eq. fl38|) is well behaved in the limit 0i — )■ 1, and thus Monte Carlo 
Simulation and thermodynamic integration may be used to calculate the value of (and 
subsequently the value of f3fi) at liquid- vapor coexistence. 
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